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Abstract 

In  this  paper,  a  two-dimensional  cross-the-channel  model  was  applied  to  investigate  the  influence  of  gas  diffusion  layer  (GDL)  property  and 
flow-field  geometry  on  the  local  reaction  rate  in  the  PEMFC  cathode  catalyst  layer.  The  model  predictions  show  that  the  rate  of  consumption 
of  oxygen  or  current  density  under  the  land  area  may  differ  considerably  from  that  under  the  channel.  Simulation  results  indicate  that  for  a 
fixed  channel  width,  increasing  the  channel-to-land  width  ratio  results  in  improved  water  transport  and  positively  impacts  the  overall  reaction 
rates  at  high  overpotential.  However,  too  high  ratio  may  retard  electron  transport  and  thus  lead  to  worse  cathode  performance.  Numerical 
simulation  results  revealed  that  GDL  electrical  conductivity  and  thickness  have  a  strong  influence  on  how  the  chemical  species  and  electrons 
are  transported  to  the  active  catalyst  sites.  Depending  on  the  electrical  conductivity  of  GDL,  the  region  of  higher  reaction  rate  may  occur  either 
under  the  land  or  under  the  channel.  Consideration  of  orthotropic  electrical  conductivity  of  GDL  affects  the  simulation  results  significantly 
highlighting  the  need  to  improve  our  understanding  of  GDL  transport  coefficients.  Simulation  of  GDL  compression  effects  showed  that  the 
total  current  density  is  not  affected  significantly  but  current  density  distribution  is.  The  extent  of  reactant  bypass  through  the  GDL  from  one 
channel  to  the  other  depends  on  the  GDL  permeability  and  results  in  a  significant  enhancement  in  reaction  rates. 

©  2004  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

PEMFCs  have  been  actively  researched  for  a  number  of 
years  because  of  their  advantages  such  as  low-temperature 
operation,  rapid  start-up  and  high  power  density  [1].  There 
is  a  continued  interest  in  improving  the  understanding  of 
various  physico-electro-chemical  processes  occurring  in  a 
PEMFC.  Of  particular  interest  are  the  cathode  processes 
which  account  for  a  majority  of  the  total  cell  electrochemical 
losses. 

The  electrochemical  reaction  occurring  at  the  cathode  of 
PEMFC  involves  three  different  species  -  ions  (protons), 
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electrons,  and  chemical  species  (oxygen  and  water).  These 
electrochemical  reactions  occur  in  a  thin  catalyst  layer,  which 
is  in  contact  with  a  proton-conducting  membrane  on  one  side 
and  a  porous  gas  diffusion  layer  (GDL)  on  the  other  side.  At 
the  cathode  side,  the  GDL  allows  the  transport  of  the  reactant 
oxygen  and  electrons  from  the  flow-field  channel  and  from 
under  the  flow-field  land,  respectively,  to  reaction  sites.  In 
addition,  the  GDL  allows  removal  of  water  produced  in  the 
catalyst  layer.  From  the  schematic  depiction  of  the  PEMFC 
cathode  in  Fig.  1,  it  can  be  readily  seen  that  the  pathways 
for  electron  and  species  transport  from  the  flow-field  plate 
to  the  catalyst  layer  are  different.  The  local  electrochemical 
reaction  rate  in  the  catalyst  layer  is,  therefore,  influenced  by 
the  geometric  parameters  and  transport  properties  of  both  the 
GDL  and  the  flow-field  plate.  Mathematical  modeling  may 
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provide  useful  insights  into  how  these  parameters  can  influ¬ 
ence  catalyst  layer  performance. 

Over  the  past  20  years,  a  wide  variety  of  PEMFC  models 
of  varying  complexity  have  been  developed.  These  include 

1- D  models  [2-7]  and  2-D  along-the-channel  (ATC)  models 
[8-14]  and,  more  recently  reported,  3-D  models  [15-19].  A 
detailed  review  of  PEMFC  models  has  been  published  by  our 
group  [20] .  A  2-D  cross-the-channel  (CTC)  model  by  West 
and  Fuller  [21]  was  applied  to  analyze  the  influence  of  the 
bipolar  plate  rib  or  land  on  local  reaction  rate.  However,  in 
their  model  the  catalyst  layer  was  treated  as  an  ultra- thin  layer, 
and  electron  transport  in  the  cathode  was  not  included.  In  fact, 
a  majority  of  the  1-D,  2-D  and  3-D  models  treat  the  catalyst 
layer  to  be  a  thin-layer  or  a  fully  flooded  thin-layer  [22-23]. 
Models  where  the  catalyst  layer  is  treated  in  the  most  detail 
are  agglomerate-type  models  [24-32],  wherein  the  catalyst 
layer  is  considered  as  a  domain  comprising  agglomerates  of 
catalyst  particles.  However,  a  majority  of  the  agglomerate 
models  are  1-D,  an  exception  being  2-D  ATC  models  [3 1-32] . 
In  summary,  owing  to  the  directions  considered  in  the  1-D  and 

2- D  ATC  models  or  because  of  the  description  of  the  catalyst 
layer  as  a  thin-film,  the  influence  of  flow-field  plate  geometry 
on  the  variation  of  the  local  electrochemical  reaction  rate  has 
not  been  reported. 

The  objective  of  this  paper  is  to  investigate  the  influence 
of  flow-field  and  GDF  geometric  parameters  and  transport 
properties  on  the  local  electrochemical  reaction  rate  distribu¬ 
tion  in  the  PEMFC  cathode  catalyst  layer.  The  catalyst  layer 
is  modeled  as  spherical-agglomerates  of  catalyst-electrolyte 
providing  detailed  distribution  of  the  reaction  rate. 


2.  Model  description 

A  detailed  description  of  the  model  employed  in  this  pa¬ 
per  is  provided  elsewhere  [33-34].  Here  we  only  briefly  de¬ 
scribe  the  computational  geometry  and  key  equations,  bound¬ 
ary  conditions,  model  assumptions,  computational  procedure 
and  input  parameters. 


Fig.  1.  Schematic  of  a  PEMFC  cathode  with  an  agglomerate  structure  cat¬ 
alyst  layer. 


2.7.  Computational  geometry 

Fig.  1  illustrates  the  cathode  structure  of  a  typical  PEMFC, 
in  which  the  catalyst  layer  is  composed  of  many  agglomer¬ 
ates.  An  agglomerate  is  made  up  of  platinum  particles,  carbon 
black  and  Nafion.  Each  agglomerate  is  covered  by  a  thin  film 
of  electrolyte.  Assuming  symmetry  in  adjacent  parallel  chan¬ 
nels  a  domain  consisting  of  a  half-channel  and  a  half  land  is 
considered  for  the  base  case.  The  rectangle  defined  by  dashed 
line  describes  the  region  of  interest.  When  reactant  bypass 
phenomenon  from  channel-to-channel  is  modeled,  the  com¬ 
putational  geometry  extends  to  v  =  xch,2-  Chemical  species 
transport  is  assumed  to  occur  by  the  sequential  pathways  - 
diffusion  through  gas  pores  in  GDF  and  catalyst  layer,  fol¬ 
lowed  by  diffusion  of  dissolved  oxygen  within  the  agglom¬ 
erate  [32]. 


2.2.  Key  Model  equations 


2.2.1.  Reactant  transport 

Chemical  species  transport  in  the  GDF  and  the  catalyst 
layer  is  described  by  the  Maxwell-Stefan  equations  [35] .  The 
diffusivities  in  the  GDF  and  catalyst  layer  were  modified 
through  the  Bruggemann  relation  [2,4]: 
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where,  ni  is  species  mass  flux  vector,  Wi  is  mass  fraction  for 
component  /,  p  is  gas  mixture  density,  M  is  molecular  weight 
of  the  gas  mixture,  Dy  is  binary  diffusivities  of  species  i  and 
j,  £  is  porosity. 

Using  Bruggeman  correlation  for  the  catalyst  layer  is  prob¬ 
ably  appropriate  because  it  was  originally  derived  for  granu¬ 
lar  packing.  However,  its  application  to  the  fibrous  GDF  may 
be  questionable.  This  topic  is  being  further  investigated  in 
our  research  group  [36]. 


2.2.2.  Electrode  reaction 

The  gaseous  oxygen  must  first  dissolve  and  then  diffuse 
through  the  electrolyte  film  surrounding  the  agglomerate, 
prior  to  participating  in  the  reaction  or  diffusing  inward  in 
the  agglomerate  while  reacting.  The  overall  reaction  rate  for 
oxygen  consumption  in  the  cathode  can  be  expressed  in  terms 
of  the  local  oxygen  concentration,  reaction  rate  constant,  an 
effectiveness  factor  for  the  electrode  reaction,  Fr,  and  ag¬ 
glomerate  physical  parameters. 

Ro2  =  El( - 1 -  +  (ragg  +  S)8  \  ~ 1  (4) 

H  \  Erkc(l  £CAT)  ^agg^agg^/ 
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where,  Ro2  is  oxygen  reaction  rate  based  on  catalyst  layer 
volume,  H  is  Henry’s  law  constant,  Ev  is  the  catalyst  effec¬ 
tiveness  factor,  ragg  is  the  agglomerate  radius,  8  is  the  thick¬ 
ness  of  the  electrolyte  film  covering  an  agglomerate,  D  is 
the  dissolved  oxygen  diffusivity  in  Nation.  In  our  model,  the 
charge-transfer  reaction  is  considered  to  be  the  rate-limiting 
step  in  the  oxygen  reduction  reaction  (ORR)  [37].  Accord¬ 
ingly,  the  electrochemical  reaction  rate  constant  (i kc )  is  given 
by  the  following  equation: 
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The  catalyst  effectiveness  factor  for  the  spherical  agglomer¬ 
ate  is  given  by  [38] 
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layer,  ionic  loss  in  the  catalyst  layer.  Thus,  the  total  potential 
loss  is  defined  as  the  nominal  cathode  overpotential  (NCO), 
which  is  the  potential  difference  between  the  electron  phase 
at  the  land  (z  =  Zg  in  Fig.  1)  and  the  electrolyte  phase  at  the  in¬ 
terface  between  catalyst  layer  and  membrane  (z  =  0  in  Fig.  1). 
This  definition  of  overpotential  holds  true  because  the  elec¬ 
tron  phase  potential  at  z  =  zg  is  considered  to  be  constant  and 
equal  to  zero.  However,  it  must  be  noted  that  the  strict  defini¬ 
tion  of  overpotential  is  the  difference  in  potential  of  the  two 
phases  at  zero  current  and  that  at  any  current. 

2.3.  Boundary  conditions 

The  computational  sub-domains  (SDs)  and  associated  di¬ 
mensions  are  shown  in  the  Fig.  2,  which  is  used  to  define 
boundary  conditions.  The  normal  fluxes  of  any  transported 
parameter,  0,  which  may  represent  species  flux,  proton  flux 
and  electron  flux,  are  assumed  to  be  zero  at  v  =  0  and  x  =  x\d, 
because  symmetry  is  assumed  at  these  boundaries. 


where  0l  is  a  dimensionless  group,  commonly  known  as 
Thiele’s  modulus  for  chemical  reactions: 


where,  De ff  is  the  effective  diffusivity  of  the  dissolved  oxygen 
in  electrolyte  phase. 

In  the  catalyst  layer,  under  electrical  load  conditions  the 
protons  are  continuously  consumed  generating  current.  Ac¬ 
cording  to  the  species  conservation,  the  divergence  of  current 
density  would  be  related  oxygen  reaction  rate: 

^  .i  =  4F^Sl( - - - h  (^agg_+3)5\ 
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(&)-n  =  0  (13) 

where,  •  n  is  normal  flux  operator. 

2.3.1.  Boundary  conditions  for  chemical  species 
transport  ( Eq.  (1)) 

Xi  =  Xito  atz  =  zg,  0  <  x  <  xch  •  (14) 

(m)  •  n  =  0  atz  =  Zg,  vch  <  x  <  x\&  and  z  =  0  :  (15) 

2.3.2.  Boundary  conditions  for  proton  charge  balance 
(Eq.  (9)) 

(iP)  •  n  =  0  atz  =  zc  '  (16) 


2.2.3.  Charge  transport 

Proton  migration  and  electron  transport  in  the  cathode  can 
be  described  by  Ohm’s  law: 

ip  =  -&fffV0i  (9) 

ie.  c  =  -kfVcps  (10) 

where,  ip  is  the  proton  flux  and  ze,c  is  the  electron  flow  in  the 
catalyst  layer,  k\  and  ks  are  proton  conductivity  and  electron 
conductivity,  respectively.  f  \  is  the  potential  in  electrolyte 
phase,  and  0S  is  the  electrical  potential  in  electron  transport. 

Conservation  of  the  charge  leads  to  the  following  equa¬ 
tions: 

in  the  catalyst  layer  :  V  •  ip  +  V  •  i  =  0  (11) 

intheGDL  :  V  •  ze,G  =  0  (12) 

where,  /e,G  is  the  electron  flow  in  the  GDL. 

Only  the  cathode  is  considered  in  this  model.  The  potential 
losses  in  the  cathode  include  activation  overpotential  at  the 
cathode,  electronic  losses  in  both  the  GDL  and  the  catalyst 
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Fig.  2.  Sub-domains  of  simulation. 
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Table  1 


Operating  and  electrode  parameters  for  base  case 


Parameters 

Value 

Units 

Sources 

Operating  pressure  (absolute) 

1.5 

atm 

Operating  temperature 

80 

°C 

Oxygen/nitrogen  ratio 

21/79 

— 

Air  feed  relative  humidity 

50% 

— 

Half  channel  width 

0.5 

mm 

Half  channel  land  width 

0.5 

mm 

GDL  thickness 

250 

|xm 

Porosity  of  GDL,  sqdl 

0.4 

— 

Pt  loading, 

0.4 

mg  cm-2 

Pt  particle  diameter,  rpt 

2.5 

nm 

Effective  specific  agglomerate  surface  area,  aa gg 

3.6  X  105 

2  -3 

m  m 

Radius  of  agglomerate,  ragg 

1 

|xm 

Catalyst  layer  thickness,  tc\ 

15 

[xm 

[28] 

Porosity  of  catalyst  layer,  scat 

0.1 

— 

[32] 

Conductivity  of  GDL,  ks 

100 

Sm-1 

[25] 

Cathodic  transfer  coefficient,  ac 

Low  slope  (>0.8  V) 

1 

— 

[50] 

High  slope  (<0.8  V) 

0.61 

— 

Reference  exchange  current  density a,  irff 

Low  slope  (>0.8  V) 

3.85  x  lO"8 

A  cm-2 

[50] 

High  slope  (<0.8  V) 

1.5  x  10“6 

A  cm-2 

O2  diffusivity  in  Nafiona,  D 

8.45  x  10“10 

2  -l 
m  s 

[50] 

Henry’s  constant,  H 

0.3125 

atm  m3  mol- 1 

[37] 

Thickness  of  electrolyte  film  covering  each  agglomerate,  8 

80 

nm 

[40] 

Electrolyte  fraction  in  agglomerate,  sa gg 

0.5 

— 

[14] 

Reference  O2  concentration13,  Cq( 

0.85 

mol  m-3 

[50] 

Effective  Pt  surface  ratio,  s\ 

0.75 

— 

[51] 

Binary  diffusivities 

PDo  2,h2o 

3.70  x  10"5 

atm  m2  s-1 

PDo2,n2 

2.79  x  10“5 

atm  m2  s-1 

[52] 

™N2,H20 

3.87  x  10“5 

atm  m2  s-1 

a  80  °C. 

b  80  °C  and  1 .5  atm  air. 


2.3.3.  Boundary  conditions  for  electron  charge  balance 
(Eq.  (10)) 

(ps  =  0  a U  =  zg,  Tch  <  x  <  x\d  :  (18) 

0e,c)  •  n  =  0  atz  =  0  :  (19) 

2.4.  Model  assumptions 

The  key  assumptions  employed  in  the  model  are  as  fol¬ 
lows: 

-  isothermal  and  steady  state  operation; 

-  partial  pressure  losses  in  the  electrode  are  only  caused  by 
multi-component  diffusion; 

-  agglomerates  are  spherical  in  shape,  mono-dispersed,  and 
evenly  distributed; 

-  local  overpotential  within  an  agglomerate  is  assumed  to  be 
constant. 

2.5.  Solution  procedure 

The  model  equations  combined  with  boundary  conditions 
listed  in  the  prior  sections  were  solved  with  the  finite  ele¬ 


ment  method  using  the  commercially  available  PDE  solver 
program  FEMLAB®.  The  stationary  non-linear  solver  op¬ 
tion,  which  is  used  for  non-linear  stationary  PDE  problems, 
was  chosen  to  solve  the  highly  non-linear  problem.  Of  three 
different  solution  forms  available,  as  recommended  in  FEM¬ 
LAB  user’s  guide  [39],  the  general  solution  form  was  chosen 
for  the  non-linear  problem. 

2.6.  Model  parameters 

The  operating  conditions  and  electrode  parameters  are 
shown  in  Table  1 . 


3.  Results  and  discussions 

Fig.  3  shows  the  cathode  polarization  curve  generated  by 
our  model  for  the  base  case.  The  cathode  potential  shown  in 
Fig.  3  is  obtained  by  subtracting  the  nominal  cathode  over¬ 
potential  (NCO)  from  the  theoretical  oxygen  reduction  po¬ 
tential.  The  NCO  includes  cathode  activation  overpotential 
electronic  losses  in  the  GDL  and  the  catalyst  layer  as  well  as 
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Fig.  3.  Cathode  polarization  curve  of  base  case. 

proton  transport  losses  in  the  electrolyte  phase  of  the  catalyst 
layer. 
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V  =  Eth  -  NCO 
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Fig.  4.  Current  density  distributions  at  z  =  0  for  three  Ch-Ld  ratios:  1:1,  2:1 
and  4:1  for  (a)  high,  (b)  intermediate  and  (c)  low  NCOs. 


Eth  =  1.229  -  8.456  x  10“4(T  -  298.15) 

+  4.31  x  10“5r  (ln[PU2]  +  hn[P02]J  (21) 

The  cathode  polarization  curve  of  Fig.  3  does  not  include 
the  potential  losses  due  to  membrane  resistance,  anode  over¬ 
potentials  and  electrical  contact  resistance.  The  polarization 
curve  follows  the  general  trend  observed  in  typical  experi¬ 
mental  results  [22,40-43]. 

3.1.  Effect  of  channel-to -land  width  ratio  on  current 
density 

One  of  the  key  flow-field  plate  design  parameters  is  the 
channel-to-land  ratio.  Whereas  a  wider  channel  may  pro¬ 
mote  chemical  species  transport,  a  corresponding  thinner 
land  width  may  restrict  electron  transport.  Because  the  elec¬ 
trochemical  reaction  involves  both  the  chemical  species  and 
electrons,  the  nature  and  extent  of  how  the  channel-to-land 
(Ch-Ld)  ratio  affects  the  cathode  performance  will  depend 
on  the  competing  effect  of  oxygen  and  electron  transport. 

To  investigate  the  effect  of  flow-field  geometric  parame¬ 
ters,  simulations  were  carried  out  for  three  Ch-Ld  ratios.  In 
all  cases,  the  channel  width  was  kept  constant  while  the  land 
width  was  varied.  Current  density  distributions  at  z  =  0  for 
three  Ch-Ld  ratios  at  NCOs  of  0.4,  0.5  and  0.6  V,  are  shown 
in  Fig.  4.  In  the  figure,  average  current  density  at  the  cata¬ 
lyst  layer/membrane  interface  (z  =  0),  /avg,  represents  the  total 


current  density  that  is  associated  with  the  net  reaction  rate  in 
the  whole  catalyst  layer. 

At  low  NCO  (case  c),  an  increase  in  Ch-Ld  ratio  does  not 
improve  the  cathode  performance.  Low  NCO  corresponds  to 
low  current  density  or  low  oxygen  consumption  rate.  Under 
these  conditions,  it  appears  that  electron  transport  in  the  cath¬ 
ode  has  a  stronger  influence  on  the  electrochemical  reaction 
rate,  so  increases  in  channel  width,  which  might  be  expected 
to  enhance  oxygen  transport,  do  not  have  any  significant  im¬ 
pact  on  cathode  performance.  At  high  NCO  (case  a),  because 
the  oxygen  consumption  rate  is  high,  species  transport  lim¬ 
itations  are  expected  to  become  important.  The  simulation 
results  show  that  at  high  NCO,  an  increase  in  Ch-Ld  ratio 
results  in  an  improvement  in  the  cathode  performance.  This 
may  be  explained  by  the  fact  that  at  high  NCO,  there  is  high 
current  density  or  high  oxygen  consumption  rate.  Thus,  oxy¬ 
gen  transport  from  the  channel  to  the  catalyst  layer  becomes 
a  dominant  factor.  By  increasing  the  channel-land  ratio,  oxy¬ 
gen  transport  to  the  region  under  the  land  is  facilitated.  For 
the  case  where  Ch-Ld  equals  2:1,  the  average  current  density 
increased  by  7%  compared  to  the  base  case.  Further  reduc¬ 
tion  in  land  width  does  not  improve  cathode  performance.  For 
the  case  of  Ch-Ld  ratio  of  4: 1,  the  current  density  is  actually 
slightly  lower  than  that  for  the  case  Ch-Ld  =  2:1.  This  may  be 
attributed  to  the  increase  in  resistance  to  electron  transport 
caused  by  supplying  electrons  to  a  larger  region  by  a  thinner 
land.  This  translates  into  reduction  in  local  overpotential  in 
the  catalyst  area  underneath  the  channel,  thereby,  a  reduction 
in  the  electrode  reaction  rate. 
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Fig.  5.  Water  vapour  partial  pressure  (atm)  distribution  in  cathode  for  three  Ch-Ld  ratios:  (a)  1:1,  (b)  2: 1  and  (c)  4:1 .  NCO  =  0.5  V. 


3.2.  Effect  of  channel- to -land  width  ratio  on  water 
transport 

A  high  Ch-Ld  ratio  is  expected  to  enhance  water  transport 
rates  and  improve  water  removal  from  catalyst  region  under 
the  channel  thereby  reducing  the  risk  of  flooding  the  catalyst 
layer.  This  phenomenon  was  investigated  by  studying  the  wa¬ 
ter  partial  pressure  distributions  in  the  cathode  for  the  three 
width  ratios.  Fig.  5  shows  water  vapor  pressure  distribution 
in  cathode  for  a  NCO  equal  to  0.5  V.  It  can  be  seen  in  the 
Fig.  5  that  the  water  partial  pressure  in  the  cathode  decreases 
significantly  with  a  reduction  in  the  land  width.  Lower  wa¬ 
ter  partial  pressure  can  be  a  result  of  either  a  less  amount 
of  water  being  produced  or  more  effective  water  removal.  In 
Fig.  4b,  the  current  density  under  the  land  area  for  Ch-Ld 
ratio  of  1 : 1  is  lower  than  that  for  Ch-Ld  ratios  of  2: 1  and  4: 1 
(and  therefore  the  water  production  rate  is  lower),  i.e.  less 
water  is  produced  under  the  land  due  to  the  cathode  reaction. 
However,  as  shown  in  Fig.  5,  the  water  vapour  pressure  under 
the  land  area  for  a  Ch-Ld  ratio  of  1 : 1  is  higher  than  that  for 
the  other  two  cases  indicating,  as  expected,  thus  a  low  Ch-Ld 
ratio  impairs  water  removal. 

It  should  be  noted  that  in  the  preceding  discussion,  any 
potential  impact  of  increased  contact  resistance  between  the 
GDL  and  the  flow-field  plates  as  the  land  width  is  reduced  has 
not  been  accounted  for.  The  impact  of  the  effect  of  the  contact 
resistance  is  currently  being  investigated  in  our  laboratory. 

3.3.  Effect  of  GDL  thickness  on  current  density 

Commercial  GDL  layers  are  available  in  variety  of  thick¬ 
nesses.  To  investigate  the  effect  of  GDL  thickness,  simula¬ 
tions  were  carried  out  for  three  GDL  thicknesses:  0.15,  0.25 
and  0.35  mm.  Fig.  6  presents  current  density  distributions  at 
the  catalyst  layer/membrane  interface  (e.g.  z  =  0)  for  three  dif¬ 
ferent  NCOs.  Average  current  densities  are  also  reported  in 
the  figure.  Intuitively,  one  might  think  that  a  thinner  GDL  will 
be  beneficial  to  the  electrode  performance,  because  it  would 
reduce  resistance  to  the  transport  of  both  oxygen  and  elec¬ 
trons.  However,  simulation  results  indicate  that  the  average 


current  density  is  not  affected  significantly  as  GDL  thickness 
is  reduced  from  0.350  to  0.15  mm. 

In  general,  increasing  the  GDL  thickness  has  the  effect  of 
flattening  the  current  density  profile  across  the  channel  and 
land  area.  However,  the  effect  of  GDL  thickness  on  local  reac¬ 
tion  rate  (current  density)  is  complicated  and  depends  on  the 
competing  effects  of  oxygen  and  electron  transport.  At  low 
overpotential,  electron  transport  plays  a  more  important  role 
in  determining  the  local  reaction  rate.  For  the  case  of  0. 15  mm 
thick  GDL,  electrons  must  travel  more  than  three  times  the 


Fig.  6.  Current  density  distributions  at  z  =  0  for  different  GDL  thicknesses: 
0.15,  0.25  and  0.35  mm  for  (a)  high,  (b)  moderate  and  (c)  low  NCOs,  re¬ 
spectively. 
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distance  from  current  collector  edge  to  a  catalyst  site  at  the 
middle  of  flow  channel  compared  to  a  catalyst  site  located 
directly  underneath  the  current  collector.  This  indicates  that 
fewer  electrons  will  participate  in  the  electrode  reaction  in 
the  region  under  the  channel.  From  the  Fig.  6c  it  can  be  seen 
that  in  the  middle  of  the  channel  area  (v  =  0),  current  den¬ 
sity  decreases  with  the  reduction  of  the  GDL  thickness.  At 
higher  current  density,  shown  in  the  Fig.  6a,  a  similar  situation 
occurs.  At  this  condition,  oxygen  diffusion  is  an  important 
rate-limiting  factor.  Similar  to  low  current  density  situation, 
the  path  for  transport  of  oxygen  under  the  land  area  is  much 
larger  than  that  under  channel  area.  As  a  result,  in  the  middle 
of  land  region  (x  =  1  mm)  the  current  density  decreases  with 
the  reduction  of  the  GDL  thickness.  In  summary,  the  cur¬ 
rent  density  distribution  depends  on  the  competitive  effect 
between  oxygen  diffusion  and  electron  transport. 

3.4.  Effect  of  GDL  conductivity 

A  wide  range  of  GDL  conductivity  values  have  been  em¬ 
ployed  in  PEMFC  models  [3,17,44-45].  To  study  the  ef¬ 
fect  of  GDL  conductivity  on  the  cathode  reaction,  simula¬ 
tions  were  carried  out  for  GDL  conductivities  of  50,  200 
and  1000 Sm'1,  in  addition  to  the  base  case  condition  of 
1 00  S  m~ 1 .  The  electron  conductivity  in  the  catalyst  layer  was 
assumed  to  be  100  S  m-1  for  all  the  cases  [53].  Fig.  7  depicts 
the  local  overpotential  distributions  in  the  catalyst  layer  for 


Fig.  8.  Current  density  distribution  at  z  =  0  for  different  GDL  conductivities. 
NCO  =  0.5  V. 

the  four  different  GDL  conductivities  at  NCO  =  0.5  V.  For 
the  case  where  the  GDL  conductivity  is  50  Sm  1 ,  the  lo¬ 
cal  overpotential  varies  up  to  20%  from  the  minimum  to  the 
maximum.  This  will  cause  a  considerable  difference  in  the 
prediction  of  local  current  density,  which  cannot  be  captured 
by  ultra-thin  catalyst  layer  models,  where  only  a  constant 
overpotential  is  used. 


u  Under  Channel 
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Fig. 7.  Local  overpotential  (V)  distributions  in  the  catalyst  layer  for  different  GDL  conductivities:  (a)  50  Sm  kOyllOOSm  1,(c)200Sm  1  and(d)1000Sm 
NCO  =  0.5  V. 
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It  should  also  be  noted  that  for  a  GDL  conductivity  of 
1000  Sm-1,  although  the  electron  transport  resistance  may 
be  negligible,  the  local  overpotentials  may  still  have  10% 
differential  locally.  This  suggests  that  proton  transport  re¬ 
sistance  in  the  catalyst  layer  also  contributes  to  lower  local 
overpotential. 


Fig.  8  illustrates  the  current  density  distributions  at  z  =  0 
for  different  GDL  conductivities  at  a  NCO  of  0.5  V.  It  can 
be  seen  that  GDL  conductivity  has  a  significant  impact  on 
current  density  prediction.  With  an  increase  in  GDL  conduc¬ 
tivity  from  50  to  lOOOSm-1,  the  predicted  current  density 
increases  by  over  30%.  This  clearly  demonstrates  the  need  of 
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Fig.  9.  Oxygen  reaction  rate  distributions  in  the  catalyst  layer  for  base  case  (isotropic  conductivity)  and  orthotropic  GDL  conductivity  for  three  NCOs:  (a) 
0.3  V,  (b)  0.5  V  and  (c)  0.7  V. 
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selecting  an  appropriate  value  for  the  GDL  conductivity  if  a 
reliable,  predictive  model  is  to  be  developed. 

In  addition,  GDL  conductivity  influences  the  current  den¬ 
sity  distribution.  The  predicted  current  density  distribution 
shown  in  Fig.  8  shows  that  with  an  increase  in  GDL  con¬ 
ductivity,  there  is  a  significant  increase  in  the  current  density 
produced  in  the  catalyst  region  underneath  the  channel  area 
and  a  modest  increase  in  the  current  density  produced  below 
the  land  area. 

3.5.  Influence  of  orthotropic  GDL  conductivity  on 
reaction  rate  distribution 

The  most  commonly  used  GDL  materials  -  carbon  pa¬ 
per  and  carbon  cloth  -  are  fibrous  by  nature  and  expected 
to  possess  conductivities  that  are  different  in  the  in-plane  di¬ 
rection  compared  to  the  through-plane  direction.  For  certain 
GDL  materials,  such  as  Toray  carbon  paper,  in-plane  con¬ 
ductivity  has  been  reported  to  be  more  than  ten  times  higher 
than  through-plane  conductivity  [46].  The  two-dimensional 
across-the-channel  model  was  used  to  investigate  the  impact 
of  the  orthotropic  variation  of  GDL  electrical  conductivity. 
For  the  simulation,  the  through-plane  conductivity  was  kept 
the  same  as  the  base  case,  and  the  in-plane  conductivity  was 
assumed  to  be  10  times  that  of  the  base  case  value.  Fig.  9 
shows  the  comparison  of  oxygen  reaction  rate  distributions 
for  the  base  case  (isotropic  conductivity)  and  orthotropic 
GDL  conductivity  case.  A  number  of  interesting  observations 
can  be  made. 

For  low  NCO  (0.3  V)  or  low  reaction  rate  conditions,  it 
can  be  observed  that  for  the  isotropic  case  the  reaction  rate 
varies  significantly  in  both  directions  -  across  the  width  and 
along  the  thickness  of  the  catalyst  layer.  The  oxygen  con¬ 
sumption  rate  is  higher  underneath  the  land  than  that  under 
the  channel.  For  the  orthotropic  case,  the  reaction  rate  varies 
primarily  along  the  thickness.  This  is  expected  because  for 
the  orthotropic  case  high  in-plane  conductivity  leads  to  a  rela¬ 
tively  uniform  overpotential  distribution  in  the  catalyst  layer. 
The  rate  of  consumption  is  low  enough  so  that  oxygen  diffu¬ 
sion  does  not  influence  the  reaction  rate  distribution  across 
the  width.  For  both  cases,  more  reaction  is  predicted  to  occur 
near  the  catalyst  layer/membrane  interface  than  that  near  the 
catalyst  layer/GDL  interface. 

Flow  Channel 


Gas 


For  intermediate  NCO  (0.5  V),  the  reaction  rate  varies  sig¬ 
nificantly  in  both  directions  -  across  the  width  and  through 
the  thickness  -  of  the  catalyst  layer  both  the  iso-  and  or¬ 
thotropic  cases.  For  the  isotropic  case,  at  any  location  along 
the  catalyst  thickness,  i.e.,  at  a  fixed  z,  the  region  of  higher 
reaction  rate  is  predicted  to  be  near  the  channel-land  inter¬ 
face.  For  the  orthotropic  case,  the  region  of  higher  reaction 
rate  is  predicted  to  be  always  underneath  the  channel  at  any 
fixed  z.  This  indicates  that  when  the  in-plane  conductivity  is 
high  oxygen  diffusion  becomes  the  limiting  process  for  the 
electrochemical  reaction.  At  both  intermediate  and  low  NCO 
conditions,  the  rate  of  oxygen  reduction  is  greater  at  the  cat¬ 
alyst  layer/membrane  interface  compared  to  at  the  catalyst 
layer/GDL  interface  for  both  the  iso-  and  orthotropic  cases. 

Interestingly,  for  a  high  NCO  (0.7  V)  the  reaction  rate 
distribution  appears  to  be  very  similar  both  the  iso-  and  or¬ 
thotropic  cases.  Further,  the  reaction  rate  is  predicted  to  vary 
primarily  across  the  channel  and  the  land,  i.e.,  the  v  direction. 
The  reaction  rate  variation  through  the  depth  of  the  catalyst 
layer,  i.e.,  the  z  direction  is  small.  Thus,  unlike  the  lower  NCO 
conditions,  the  rate  of  oxygen  consumption  is  no  longer  con¬ 
siderably  higher  near  the  catalyst  layer/membrane  interface. 

3.6.  GDL  compression  effect 

In  an  assembled  fuel  cell  stack,  the  GDL  thickness  varies 
depending  on  the  location.  Usually,  the  portion  under  land 
will  be  compressed  by  approximately  15%  [47].  This  com¬ 
pression  will  probably  affect  the  electronic  conductivity  and 
the  chemical  species  transport.  A  simulation  was  conducted 
to  investigate  how  local  variation  of  conductivity  and  poros¬ 
ity  in  the  GDL  caused  by  compression  might  affect  the  local 
oxygen  consumption  rate  and,  hence,  the  current  density. 

For  the  simulation,  the  GDL  under  land  area  was  as¬ 
sumed  to  be  15%  thinner  than  that  under  the  channel  area 
(see  Fig.  10).  As  a  result  of  compression,  the  porosity  and 
electronic  conductivity  of  the  GDL  under  the  land  region 
were  assumed  to  decrease  and  increase  by  15%,  respectively 
with  respect  to  the  base  case  values.  The  change  in  binary 
diffusivities  due  to  porosity  change  was  accounted  for  by  Eq. 
(2). 

Fig.  1 1  illustrates  the  current  density  distribution  at  z  =  0 
for  the  base  case  and  with  GDL  compression.  It  can  be  seen 
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Fig.  10.  Cathode  structure  with  compressed  GDL. 
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Fig.  11.  Current  density  distributions  at  z  =  0  for  base  case  and  GDL  com¬ 
pression  case.  NCO  =  0.5  V. 

that  for  a  15%  increase  in  conductivity  and  a  15%  decrease  in 
porosity  the  overall  cathode  performance  is  not  affected  sig¬ 
nificantly.  The  fact  that  the  two  simultaneous  changes  appear 
to  have  a  very  small  net  effect  may  be  because  the  oxygen 
reaction  rate  is  determined  by  two  key  factors  -  local  oxygen 
concentration  and  local  overpotential.  A  compressed  GDL 
increases  the  resistance  to  oxygen  diffusion  to  the  region  un¬ 
der  land,  resulting  in  a  lower  local  oxygen  concentration. 
However,  at  the  same  time  compression  increases  the  local 


conductivity  and  it  becomes  easier  to  transport  electrons  to 
the  region  under  the  land,  thus  increasing  local  overpotential. 
Therefore,  one  effect  appears  to  be  balanced  by  the  other. 
The  compression  effect,  however,  may  become  more  impor¬ 
tant  in  a  situation  where  a  stack  is  over-compressed  or  there 
is  variability  in  humidification  resulting  in  the  GDL  being 
deformed  with  a  concentration  loss  in  species  transport  and 
no  compensatory  increase  in  electrical  conductivity.  Further 
investigation  of  these  effects  is  planned. 

Although  the  predicted  total  current  density  does  not  vary 
significantly  with  the  GDL  compression,  the  current  den¬ 
sity  distribution  changes  noticeably.  Fig.  1 1  shows  that  more 
current  is  produced  under  the  channel  area,  because  of  the 
reduced  oxygen  diffusion  in  compressed  region. 

3. 7.  Reactant  bypass  phenomena 

In  a  fuel  cell  flow-field  plate  with  serpentine  type  chan¬ 
nels,  the  reactant  partial  pressure  as  well  as  total  pressure 
decreases  along  the  flow  path.  Accordingly,  there  are  certain 
locations  in  the  flow-field  plate  where  both  the  reactant  partial 
pressure  and  total  pressure  in  the  adjacent  flow-channels  dif¬ 
fer  significantly.  These  conditions  promote  diffusive  and/or 
convective  transport  of  reactants  through  the  GDL  from  one 
channel  to  a  neighboring  channel  where  partial  and  total  pres¬ 
sure  is  lower.  This  phenomenon  is  often  referred  as  reactant 


Nominal  Cathode  Overpotential  =G.25V 


Nominal  Cahtode  Overpotential  =  0.5V 


Nominal  Cathode  Overpotential=0.25V 


Nominal  Cathode  Overpotential=0.5V 


Channel(L) 


Current  Collector 


;hannel(R) 


I 


(C) 


0 


‘  \ 
»  \ 
l  \ 
1  \ 


A  ^ 
A  A 
A  A 


A  A 
A  A 
A  A 
A  A 
A  A 
A  A 
A  A 
A  A 
A  A 
A  A 
A  A 


A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 


A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 


0.5 


1.5 


x  10 


2 

■3 


Channel(L) 


(d) 


o 


. 


Current  Collector 


A  •* 
A  s 

\  ^ 

V  ^ 

h 

V  \ 

11 


\ 

\ 

i 


i  t 


-  e 
'  / 
'  / 
'  / 


'  \ 
I 

i 


Channel(R) 


LI 


0.5 


1.5 


x  10 


2 

■3 


Fig.  12.  Oxygen  flux  in  the  GDL.  High  permeability  cases  under  NCOs:  (a)  0.25  V  and  (b)  0.5  V.  Low  permeability  cases  under  NCOs:  (c)  0.25  V  and  (d)  0.5  V. 
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Fig.  13.  The  effect  of  convective  flow  between  adjacent  channels  on  current  density  distribution  at  the  electrode-membrane  interface.  NCOs:  (a)  0.25  V  and 
(b)  0.5  V. 


bypass  [48].  Mench  et  al.  [49]  have  suggested  that  this  “unde¬ 
sired”  reactant  penetration  may  lead  to  dramatically  reduced 
performance  for  fuel  cell  due  to  the  reactant  starvation.  The 
two-dimensional  cross-the-channel  model  was  used  to  inves¬ 
tigate  the  impact  of  reactant  bypass  on  cathode  performance. 

The  convective  flow  in  the  porous  GDL  between  adjacent 
section  of  channel  may  be  described  by  Darcy’s  law: 

k^)  VP  (22) 

\1  ) 

where,  u  is  the  velocity  of  the  convective  flow,  kverm  is  the 
gas  permeability  in  the  porous  media,  /x  is  the  gas  dynamic 
viscosity,  P  is  the  local  gas  pressure. 

The  continuity  equation  can  be  described  by  the  following 
equation: 

V  •  ( pu )  =  0  (23) 

where,  p  is  the  local  gas  density. 

Simulations  were  conducted  for  conditions  where  the  total 
pressure  and  oxygen  partial  pressure  in  an  adjacent  channel 
were  arbitrarily  assigned  to  be  90%  of  those  in  the  base  case 
channel.  This  assumption  may  not  necessarily  represent  a 
real  condition;  but,  it  allows  us  to  investigate  the  importance 
of  convective  flow  in  the  GDL.  For  all  the  simulations,  the 
dynamic  viscosity  of  gas,  /z,  was  kept  constant  at  2  x  10-5 
[Pa  s] .  While,  because  a  wide  range  of  gas  permeability  in 
GDL,  kperm,  ranging  from  10-11  to  10-19m2  has  been  re¬ 
ported  in  the  literature,  two  permeability  values,  2.3  x  10“ 11 
[18]  and  10“ 15  [m2]  were  used  to  account  for  the  effect  of 
the  parameter.  In  order  to  investigate  the  influence  of  the  con¬ 
vective  flow  generated  by  total  pressure  difference  between 
the  neighboring  channels  on  the  cathode  performance,  we 
also  simulated  the  case  without  the  consideration  of  the  con¬ 
vective  flow,  where  only  oxygen  partial  pressure  difference 
between  two  channels  was  considered  but  the  total  pressure 
was  assumed  to  be  constant.  This  case  can  be  treated  as  a 
special  case  where  the  gas  permeability  is  equal  to  zero. 


Fig.  12  shows  the  oxygen  flux  in  the  GDL  for  two  different 
NCOs.  When  high  gas  permeability  is  used,  as  shown  in  the 
cases  (a)  and  (b),  there  is  significant  oxygen  flux  from  one 
channel  to  the  other  and  towards  the  catalyst  layer.  This  leads 
to  an  increase  in  oxygen  concentration  underneath  the  land 
area.  For  low  gas  permeability  in  the  GDL,  when  the  NCO  is 
low  (see  Fig.  12c),  a  portion  of  the  oxygen  flux  is  also  directed 
towards  the  adjacent  channel.  When  the  NCO  increases,  the 
driving  force  for  oxygen  transport  to  the  catalyst  layer  may 
be  higher  than  that  for  the  transport  of  oxygen  to  the  adjacent 
channel  due  to  the  total  and  partial  pressure  difference  of 
oxygen.  As  such,  there  is  no  net  oxygen  flux  between  adjacent 
channels  (see  Fig.  12d). 

To  illustrate  the  effect  of  convective  flow  under  the  land, 
the  predicted  current  density  at  z  =  0  for  both  with  and  without 
convective  flow  are  shown  in  Fig.  13.  There  is  a  significant 
enhancement  in  the  electrochemical  reaction  under  the  land 
area  for  high  permeability  cases,  indicating  that  the  reactant 
bypass  phenomenon  is  beneficial  to  cathode  performance.  For 
low  permeability  conditions  the  convective  flow  has  smaller 
impact  on  the  cathode  performance,  but  the  predicted  current 
density  is  not  reduced  due  to  reactant  bypass.  Overall,  the 
simulation  results  indicate  that  reactant  bypass  will  have  a 
positive  effect  on  cathode  performance  contrary  to  previous 
published  predictions  [49] . 

4.  Conclusions 

The  influence  of  geometric  parameters  and  transport  co¬ 
efficients  of  the  gas  diffusion  layer  and  flow-field  plate  on 
the  cathode  performance  has  been  investigated  using  an 
improved  two-dimensional  electrode  catalyst  agglomerate 
model. 

Simulations  done  to  study  the  effect  of  channel-to-land 
width  ratio  indicate  that  increasing  Ch-Ld  ratio  results  in  im¬ 
proved  water  transport  and  positively  affects  the  overall  oxy¬ 
gen  reduction  reaction  rate  at  high  overpotential.  Simulation 
results  using  an  improved  model  of  the  GDL  indicate  that 
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the  transport  for  both  gas  and  electrons  plays  an  important 
role  in  cathode  performance.  Both  the  electronic  conductivity 
and  GDL  thickness  could  be  key  parameters  to  be  optimized 
to  facilitate  transport  of  electrons  and  chemical  species.  Dif¬ 
ferences  between  in-  and  through-plane  conductivities  in  the 
GDL  had  a  significant  impact  the  reaction  rate  distribution  in 
the  catalyst  layer. 

The  improved  two-dimensional  cross-the-channel  model 
was  further  used  to  conduct  preliminary  studies  of  the  ef¬ 
fect  of  GDL  compression  and  channel-to-channel  reactant 
bypass.  The  simulation  results  indicate  that  moderate  GDL 
compression  does  not  significantly  influence  cathode  perfor¬ 
mance  for  single-phase  flow  conditions.  Moreover,  channel- 
to-channel  bypass  flow  under  the  land  is  predicted  to  have  a 
positive  impact  on  cathode  performance. 
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